Perturbation Trasmission and Graph visualization

# Import necessary libraries
import networkx as nx  # For creating and manipulating graphs
import numpy as np  # For numerical operations
import matplotlib.pyplot as plt  # For plotting
from matplotlib.patches import Patch  # For creating custom legend patches

# Set the plotting style to 'ggplot'
plt.style.use('ggplot')

# Set the default figure size to [10, 10]
plt.rcParams['figure.figsize'] = [10, 10]

SIRS model simulation

def dW(delta_t: float) -> float:
    """Sample a random number at each call."""
    return np.random.normal(loc=0.0, scale=np.sqrt(delta_t))

def Euler_Maruyama_method(parameters):
    """
    Perform Euler-Maruyama method to simulate an epidemiological model.

    Parameters:
    - parameters (dict): Dictionary containing model parameters.
        Expected keys:
        't_in': Start time of simulation.
        't_end': End time of simulation.
        'N': Number of time steps.
        'mu': Natural birth and death rate.
        'b0': Transmission parameter.
        'b1': Transmission parameter.
        'phi': Transmission parameter.
        'gamma': Recovery rate.
        'ni': Disease-induced death rate.
        'alpha': Scaling factor for stochastic term.
        'S_in': Initial susceptible population.
        'I_in': Initial infectious population.
        'R_in': Initial recovered population.

    Returns:
    - TS (numpy.ndarray): Array of time steps.
    - Ss (numpy.ndarray): Array of susceptible population over time.
    - Is (numpy.ndarray): Array of infectious population over time.
    - Rs (numpy.ndarray): Array of recovered population over time.
    """
    t_in = parameters['t_in']
    t_end = parameters['t_end']
    N = parameters['N']
    mu = parameters['mu']
    b0 = parameters['b0']
    b1 = parameters['b1']
    phi = parameters['phi']
    gamma = parameters['gamma']
    ni = parameters['ni']
    alpha = parameters['alpha']
    S_in = parameters['S_in']
    I_in = parameters['I_in']
    R_in = parameters['R_in']

    # Calculate time step size
    dt = float((t_end - t_in) / N)
    # Create array of time steps
    TS = np.arange(t_in, t_end + dt, dt)
    assert TS.size == N + 1

    # Initialize arrays for S, I, R populations
    Ss = np.zeros(TS.size)
    Is = np.zeros(TS.size)
    Rs = np.zeros(TS.size)

    # Set initial values
    Ss[0] = S_in
    Is[0] = I_in
    Rs[0] = R_in

    # Perform Euler-Maruyama method to simulate SIR model
    for i in range(1, TS.size):
        t = t_in + (i - 1) * dt
        S = Ss[i - 1]
        I = Is[i - 1]
        R = Rs[i - 1]

        # Calculate effective transmission rate with stochastic term
        b0_tilde = b0 + alpha * dW(dt)
        beta = b0_tilde * (1 + b1 * np.cos(2 * np.pi * t + phi))

        # Calculate stochastic term
        dW_dt = dW(dt)

        # Update S, I, R populations using differential equations
        Ss[i] = S + ((mu - mu * S - beta * S * I + gamma * R) * dt - (alpha / b0_tilde) * beta * S * I * dW_dt)
        Is[i] = I + ((beta * S * I - ni * I - mu * I) * dt - (alpha / b0_tilde) * beta * S * I * dW_dt)
        Rs[i] = R + (ni * I - mu * R - gamma * R) * dt

    return TS, Ss, Is, Rs

# Example parameters dictionary
parameters = {
    't_in': 0,
    't_end': 1,
    'N': 5000,
    'mu': 0.009,
    'b0': 36.4,
    'b1': 0.38,
    'phi': 1.07,
    'gamma': 1.8,
    'ni': 36,
    'alpha': 0.25,
    'S_in': 0.998,
    'I_in': 0.002,
    'R_in': 0.0
}

# Run simulation using Euler-Maruyama method
TS, Ss, Is, Rs = Euler_Maruyama_method(parameters)

# Calculate basic reproduction number R0
def calculate_R0(parameters):
    """
    Calculate the basic reproduction number R0 for an SIR-type epidemiological model.

    Parameters:
    - parameters (dict): Dictionary containing model parameters.
        Expected keys:
        'mu': Natural birth and death rate.
        'gamma': Recovery rate.
        'ni': Disease-induced death rate.
        'b0': Transmission parameter.
        'phi': Transmission parameter.
        'b1': Transmission parameter.

    Returns:
    - R0 (float): Basic reproduction number.
    """
    mu = parameters['mu']
    gamma = parameters['gamma']
    ni = parameters['ni']
    beta = parameters['b0'] * parameters['phi'] * parameters['b1']
    R0 = beta / (mu + ni + gamma)
    return R0

# Calculate R0
R0 = calculate_R0(parameters)

# Print basic reproduction number (R0)
print(f"Basic Reproduction Number (R0): {round(R0, 1)}")

# Interpretation of R0
if round(R0) < 1:
    print("The infection will likely die out over time.")
elif round(R0) == 1:
    print("The infection will remain stable in the population but not cause an outbreak.")
else:
    print("The infection will likely spread and cause an epidemic.")

# Additional statistics
peak_infection = max(Is)
final_size = Ss[-1] + Is[-1] + Rs[-1]
print(f"Peak Infection: {peak_infection}")
print(f"Final Size of the Epidemic: {final_size}")

Basic Reproduction Number (R0): 0.4
The infection will likely die out over time.
Peak Infection: 0.0030332320384369007
Final Size of the Epidemic: 1.0009432020818723

SIRS PLOTS

# Function to create a plot
def plot_population(TS, data, label, color):
    plt.figure(figsize=(10, 6))
    plt.plot(TS, data, label=label, color=color)
    plt.xlabel('Time')
    plt.ylabel('Proportion')
    plt.title(f'{label} Population Over Time')
    plt.legend()
    plt.grid(True)
    plt.show()

# Plot Susceptible over time
plot_population(TS, Ss, 'Susceptible', 'blue')

# Plot Infected over time
plot_population(TS, Is, 'Infected', 'red')

# Plot Recovered over time
plot_population(TS, Rs, 'Recovered', 'green')

# Plot all three (S, I, R) on the same graph
plt.figure(figsize=(10, 6))
plt.plot(TS, Ss, label='Susceptible', color='blue')
plt.plot(TS, Is, label='Infected', color='red')
plt.plot(TS, Rs, label='Recovered', color='green')
plt.xlabel('Time')
plt.ylabel('Proportion')
plt.title('SIR Model Over Time')
plt.legend()
plt.grid(True)
plt.show()

NetworkX

1. Small-World Networks (Watts-Strogatz Model)

Small-world networks have high clustering like regular networks but also have short average path lengths between nodes. The Watts-Strogatz model can be used to generate small-world networks.

np.random.seed(42)

# Generate a Watts-Strogatz small-world network
n = 300
k = 4  # Each node is connected to k nearest neighbors
p = 0.15  # Probability of rewiring each edge
G = nx.watts_strogatz_graph(n, k, p)

# Ensure there is at least one initially infected node
initial_states = np.random.choice(['S', 'I', 'R'], size=len(G.nodes), p=[parameters['S_in'], parameters['I_in'], parameters['R_in']])
while not any(state == 'I' for state in initial_states):
    initial_states = np.random.choice(['S', 'I', 'R'], size=len(G.nodes), p=[parameters['S_in'], parameters['I_in'], parameters['R_in']])
nx.set_node_attributes(G, dict(zip(G.nodes, initial_states)), 'state')

# Define colors for each state and create a custom legend
state_colors = {'S': 'blue', 'I': 'red', 'R': 'green'}
legend_elements = [Patch(facecolor=color, edgecolor=color, label=state) for state, color in state_colors.items()]

# Fixed layout for consistent node positions in plots
pos = nx.spring_layout(G)

# Function to update node states based on SIRS dynamics
def update_states(G, t, Is, Rs, Ss):
    for node in G.nodes:
        state = G.nodes[node]['state']
        if state == 'S' and np.random.random() < Is[t]:
            G.nodes[node]['state'] = 'I'
        elif state == 'I' and np.random.random() < Rs[t]:
            G.nodes[node]['state'] = 'R'
        elif state == 'R' and np.random.random() < Ss[t]:
            G.nodes[node]['state'] = 'S'

# Simulation and plotting
infection_alive = True  # Flag to track if infection is still spreading

# Euler-Maruyama method simulation
TS, Ss, Is, Rs = Euler_Maruyama_method(parameters)

# Determine plot intervals
plot_intervals = np.linspace(0, len(TS) - 1, 10).astype(int)

# Simulation and plotting
for t_index in plot_intervals:
    # Update node colors based on their states
    node_color = [state_colors[G.nodes[node]['state']] for node in G.nodes]

    # Draw the graph with updated node colors
    plt.figure(figsize=(10, 6))
    nx.draw(G, pos, with_labels=False, node_color=node_color, node_size=50)
    plt.legend(handles=legend_elements, loc='best')
    plt.title(f'Small-World Network (Watts-Strogatz Model) at t={TS[t_index]:.2f}')
    plt.show()

    # Determine the range for the next state update
    next_t_index = plot_intervals[np.searchsorted(plot_intervals, t_index) + 1] if t_index != plot_intervals[-1] else len(TS)
    
    # Update states of nodes in the specified range
    for t in range(t_index, next_t_index):
        update_states(G, t, Is, Rs, Ss)
        
        # Check if all nodes are in susceptible state ('S')
        if all(G.nodes[node]['state'] == 'S' for node in G.nodes):
            infection_alive = False

    # If infection has died out, print message and break out of loop
    if not infection_alive:
        print("The infection has died out.")
        break

# Calculate effective reproduction number Rt over time
Rt = calculate_R0(parameters) * Ss / len(G.nodes)
print(f"Effective Reproduction Number (Rt) at the beginning: {Rt[0]}")

The infection has died out.
Effective Reproduction Number (Rt) at the beginning: 0.0013022154795595405

Various metrics to analyze the graph's properties

# Degree distribution
degree_sequence = [d for n, d in G.degree()]

# Average degree
avg_degree = sum(degree_sequence) / len(degree_sequence)

print(f"Average degree: {avg_degree}")

# Clustering coefficient
avg_clustering = nx.average_clustering(G)

print(f"Average clustering coefficient: {avg_clustering}")

# Average shortest path length
avg_shortest_path_length = nx.average_shortest_path_length(G)

print(f"Average shortest path length: {avg_shortest_path_length}")

# Diameter (may need to handle exceptions if graph is not connected)
try:
    diameter = nx.diameter(G)
    print(f"Diameter: {diameter}")
except nx.NetworkXError:
    print("Graph is not connected.")

# Degree centrality
degree_centrality = nx.degree_centrality(G)
# Example: Print degree centrality of the first few nodes
print("\nDegree Centrality\n")
for node in list(G.nodes())[:5]:
    print(f"Node {node}: Degree Centrality = {degree_centrality[node]}")

# Betweenness centrality
betweenness_centrality = nx.betweenness_centrality(G)

# Example: Print betweenness centrality of the first few nodes
print("\nBetweenness Centrality\n")
for node in list(G.nodes())[:5]:
    print(f"Node {node}: Betweenness Centrality = {betweenness_centrality[node]}")

# Assortativity
assortativity = nx.degree_assortativity_coefficient(G)
print(f"Degree Assortativity: {assortativity}")


# Global efficiency
global_efficiency = nx.global_efficiency(G)
print(f"Global Efficiency: {global_efficiency}")


# Local efficiency
local_efficiency = nx.local_efficiency(G)
print(f"Local Efficiency: {local_efficiency}")

# Degree histogram
degree_sequence = [d for n, d in G.degree()]
plt.hist(degree_sequence, bins='auto', alpha=0.75)
plt.xlabel('Degree')
plt.ylabel('Frequency')
plt.title('Degree Distribution')
plt.show()

# Community detection using Louvain method
import community

partition = community.best_partition(G)
print(f"Number of communities detected: {len(set(partition.values()))}")

# Community detection using Louvain method
partition = community.best_partition(G)
print(f"Number of communities detected: {len(set(partition.values()))}")

# Modularity
modularity = community.modularity(partition, G)
print(f"Modularity: {modularity}")

# Visualization of communities (optional)
pos = nx.spring_layout(G)
plt.figure(figsize=(12, 12))
nx.draw_networkx_nodes(G, pos, node_size=20, cmap=plt.cm.RdYlBu, node_color=list(partition.values()))
nx.draw_networkx_edges(G, pos, alpha=0.4)
plt.title('Graph with Communities (Louvain method)')
plt.show()

Average degree: 4.0
Average clustering coefficient: 0.31297142857142857
Average shortest path length: 5.746570281124498
Diameter: 12

Degree Centrality

Node 0: Degree Centrality = 0.020080321285140562
Node 1: Degree Centrality = 0.01606425702811245
Node 2: Degree Centrality = 0.020080321285140562
Node 3: Degree Centrality = 0.028112449799196783
Node 4: Degree Centrality = 0.020080321285140562

Betweenness Centrality

Node 0: Betweenness Centrality = 0.02705176631312766
Node 1: Betweenness Centrality = 0.017062091061331035
Node 2: Betweenness Centrality = 0.028341754057983386
Node 3: Betweenness Centrality = 0.07175942440050946
Node 4: Betweenness Centrality = 0.01647276093760317
Degree Assortativity: 0.029907628017542705
Global Efficiency: 0.20521358507622173
Local Efficiency: 0.3944952380952381

Number of communities detected: 16
Number of communities detected: 15
Modularity: 0.745494

2. Random graph (Erdős-Rényi)

np.random.seed(42)

# Create an Erdős-Rényi random graph
G = nx.erdos_renyi_graph(n=200, p=0.15)  # n: number of nodes, p: probability of edge creation

# Assign initial states to nodes based on S_in, I_in, R_in
initial_states = np.random.choice(['S', 'I', 'R'], size=len(G.nodes), p=[parameters['S_in'], parameters['I_in'], parameters['R_in']])
for node, state in zip(G.nodes, initial_states):
    G.nodes[node]['state'] = state

# Simulate SIRS dynamics on the graph
TS, Ss, Is, Rs = Euler_Maruyama_method(parameters)

pos = nx.spring_layout(G)  # Fixed layout for all iterations

# Determine the indices of the time points to plot
plot_intervals = np.linspace(0, len(TS) - 1, 10).astype(int)

# Define colors for each state
node_colors = {'S': 'blue', 'I': 'red', 'R': 'green'}

# Create custom legend
legend_elements = [
    Patch(facecolor='blue', edgecolor='blue', label='Susceptible (S)'),
    Patch(facecolor='red', edgecolor='red', label='Infected (I)'),
    Patch(facecolor='green', edgecolor='green', label='Recovered (R)')
]

for t_index in plot_intervals:
    node_colors = {'S': 'blue', 'I': 'red', 'R': 'green'}
    node_color = [node_colors[G.nodes[node]['state']] for node in G.nodes]

    # Draw the graph with updated node colors
    nx.draw(G, pos, with_labels=True, node_color=node_color)
    plt.legend(handles=legend_elements, loc='best')
    plt.title(f'Random graph (Erdős-Rényi) at t={TS[t_index]:.2f}')
    plt.show()

    

    # Update states of nodes based on SIRS simulation results up to the next plot interval
    for t in range(t_index, len(TS) if t_index == plot_intervals[-1] else plot_intervals[plot_intervals.tolist().index(t_index) + 1]):
        for node in G.nodes:
            if G.nodes[node]['state'] == 'S':
                if np.random.random() < Is[t]:  # Infection dynamics
                    G.nodes[node]['state'] = 'I'
            elif G.nodes[node]['state'] == 'I':
                if np.random.random() < Rs[t]:  # Recovery dynamics
                    G.nodes[node]['state'] = 'R'
            elif G.nodes[node]['state'] == 'R':
                if np.random.random() < Ss[t]:  # Loss of immunity dynamics
                    G.nodes[node]['state'] = 'S'

# Calculate effective reproduction number Rt over time
Rt = R0 * Ss / len(G.nodes)
print(f"Effective Reproduction Number (Rt) at the beginning: {Rt[0]}")

Effective Reproduction Number (Rt) at the beginning: 0.005045444444444444

3. Scale-Free Networks (Barabási-Albert Model)

Scale-free networks are characterized by a small number of highly connected nodes (hubs) and many nodes with only a few connections. The Barabási-Albert model is a popular method to generate scale-free networks.

np.random.seed(42)

# Generate a Barabási-Albert scale-free network
n = 200
m = 3  # Number of edges to attach from a new node to existing nodes
G = nx.barabasi_albert_graph(n, m)

# Assign initial states to nodes based on S_in, I_in, R_in
initial_states = np.random.choice(['S', 'I', 'R'], size=len(G.nodes), p=[parameters['S_in'], parameters['I_in'], parameters['R_in']])
for node, state in zip(G.nodes, initial_states):
    G.nodes[node]['state'] = state

# Simulate SIRS dynamics on the graph
TS, Ss, Is, Rs = Euler_Maruyama_method(parameters)

# Simulate SIRS dynamics on the graph
pos = nx.spring_layout(G)  # Fixed layout for all iterations

# Determine the indices of the time points to plot
plot_intervals = np.linspace(0, len(TS) - 1, 10).astype(int)

# Define colors for each state
node_colors = {'S': 'blue', 'I': 'red', 'R': 'green'}

# Create custom legend
legend_elements = [
    Patch(facecolor='blue', edgecolor='blue', label='Susceptible (S)'),
    Patch(facecolor='red', edgecolor='red', label='Infected (I)'),
    Patch(facecolor='green', edgecolor='green', label='Recovered (R)')
]

for t_index in plot_intervals:
    node_colors = {'S': 'blue', 'I': 'red', 'R': 'green'}
    node_color = [node_colors[G.nodes[node]['state']] for node in G.nodes]

    # Draw the graph with updated node colors
    nx.draw(G, pos, with_labels=True, node_color=node_color)
    plt.legend(handles=legend_elements, loc='best')
    plt.title(f'Scale-Free Networks (Barabási-Albert Model) at t={TS[t_index]:.2f}')
    plt.show()

    

    # Update states of nodes based on SIRS simulation results up to the next plot interval
    for t in range(t_index, len(TS) if t_index == plot_intervals[-1] else plot_intervals[plot_intervals.tolist().index(t_index) + 1]):
        for node in G.nodes:
            if G.nodes[node]['state'] == 'S':
                if np.random.random() < Is[t]:  # Infection dynamics
                    G.nodes[node]['state'] = 'I'
            elif G.nodes[node]['state'] == 'I':
                if np.random.random() < Rs[t]:  # Recovery dynamics
                    G.nodes[node]['state'] = 'R'
            elif G.nodes[node]['state'] == 'R':
                if np.random.random() < Ss[t]:  # Loss of immunity dynamics
                    G.nodes[node]['state'] = 'S'


# Calculate effective reproduction number Rt over time
Rt = R0 * Ss / len(G.nodes)
print(f"Effective Reproduction Number (Rt) at the beginning: {Rt[0]}")

Effective Reproduction Number (Rt) at the beginning: 0.005045444444444444